/*
* Tree.cpp
*
*  Created on: Mar 28, 2011
*      Author: hxu
*/

#include "Btree.h"
#include<stdlib.h>
#include<limits>
#include<cmath>
#include "Toolset.h"
#include<map>
#include "TSV.h"
#include "Wire.h"

using namespace Cxuh;
using namespace std;

namespace Cxuh {


	int Btree::TSV_bound_per_tier = 0;
	int Btree::tier_of_src = 0;
	const double Btree::ELMORE_DISTR = 0.38; // distributed Elmore coefficients
	const double Btree::ELMORE_LUMP = 0.69; // lumped Elmore coefficients

	Btree::Btree(int idd) :
	Network<int>() {
		// Auto-generated constructor stub
		root = NULL;
		id = idd;
		this->nsink = 0;
		num_mp = 0;
		num_tsv = 0;
		srcz = 1; // srcz is assumed to be put on the first tier in the beginning
		num_buf = 0;
		pd_min = numeric_limits<double>::max();
		pd_max = numeric_limits<double>::min();
	}

	Btree::~Btree() {
		// Auto-generated destructor stub
		deleteNode(root);
	}

	void Btree::deleteNode(BtreeNode* nd) {
		// free memory
		if (!nd) {
			deleteNode(nd->lchild);
			deleteNode(nd->rchild);
			delete nd;
		}
	}

	void Btree::addSink(BtreeNode* newsink) {
		// add an element to the end of the children of the root
		sinks.push_back(newsink);
		nsink = nsink + 1;
	}

	ostream& operator <<(ostream& out, Btree& tree1) {
		// output the content of trees
		BtreeNode* root = tree1.getRoot();
		int ns = tree1.nsink;
		int xmin, xmax, ymin, ymax;
		double dmin, dmax;

		out << "NumSinks: " << ns << endl;
		out << "NumBufs: " << tree1.num_buf << endl;
		out << "NumTSVs: " << tree1.num_tsv << endl;
		out << "NumMPs: " << tree1.num_mp << endl;
		tree1.counter = 0;

		/* output tree nodes */
		if (root == NULL)
			out << "Not synthesized yet." << endl;
		else {
			tree1.findXYMinMax(xmin, xmax, ymin, ymax);
			out << endl;
			out << "X: " << xmin << " " << xmax << endl;
			out << "Y: " << ymin << " " << ymax << endl;
			out << "Source : " << tree1.getRoot()->getX() << " "
				<< tree1.getRoot()->getY() << " " << tree1.getRoot()->getZ()
				<< endl << endl;
			tree1.traverse(out, root);
			out << "# Finish traversing the tree." << endl;
		}

		// output to the screen
		if (ns != tree1.counter)
			cerr << "Nsink = " << ns << " but counter = " << tree1.counter << endl;

		dmin = tree1.get_pathd_min();
		dmax = tree1.get_pathd_max();

		cout << " ---------------------------------" << endl;
		cout << " * Summary" << endl;
		cout << "   sinks: " << ns << endl;
		cout << "   merge points: " << tree1.getNumMP() << endl;
		cout << "   TSVs: " << tree1.getNumTSV() << endl;
		cout << "   buffers: " << tree1.getNumBuf() << endl;
		cout << "   Min delay [s]: " << dmin << endl;
		out << "# Min delay [s]: " << dmin << endl;
		cout << "   Max delay [s]: " << dmax << endl;
		out << "# Max delay [s]: " << dmax << endl;
		cout << "   Skew [s] : " << dmax - dmin << endl;
		out << "# Skew [s]: " << dmax - dmin << endl;
		return out;
	}

	void Btree::traverse(ostream& out, BtreeNode* rt) {
		// preorder traverse
		BtreeNode* lchild = rt->lchild;
		BtreeNode* rchild = rt->rchild;
		BtreeNode* parent = rt->parent;

		if (rt->ntype == LEAF)
			counter += 1;

		// output preceding nets first
		if (parent) {
			int znet = rt->getZ(); // net is in the same tier with its sink
			int x0, x1, y0, y1;
			x0 = parent->getX();
			x1 = rt->getX();
			y0 = parent->getY();
			y1 = rt->getY();

			// zero-length net is also recorded
			out << "Net : " << parent->name << "_" << rt->name << endl; // net id equals to Node id, may be duplicated
			out << "\tCoordinate : " << x0 << " " << y0 << " " << x1 << " " << y1
				<< " " << znet << endl;

		}

		// output rt
		out << *rt;

		// output children
		if (lchild)
			traverse(out, lchild);
		if (rchild)
			traverse(out, rchild);
	}

	int Btree::getID() const {
		return id;
	}

	void Btree::CTS3D() {
		// clock tree synthesis, top-down with MMM3D, bottom-up with DME3D
		srcz = tier_of_src;
		root = MMM3D(sinks, TSV_bound_per_tier, 'x'); // generate topology
		root->z = srcz; // change to source z
		//DME3D();
		DME3D_simple(root, true); // build the tree w or w/o buffers
		//buffering(root); // separately insert buffers bottom up
		update();	// update the count and label of all nodes
	}

	void Btree::DME3D_simple(BtreeNode *rt, const bool insert_buf){
		// simplified version of DME3D, without merging segments
		if (!rt)
			throw(" Error: DME3D_simple(NULL)!");

		double Croot;
		double dmax;
		dmax = DME3D_simple(rt, Croot, insert_buf);
		cout << " * Before buffering, largest delay is : " << dmax << " s."<< endl;
	}

	double Btree::DME3D_simple(BtreeNode *rt, double& c, const bool insert_buf){
		// sub procedure of DME3D_simple(rt), get the load C and the largest delay from rt to the sinks
		if (!rt)
			throw(" Error: DME3D_simple(NULL, c)!");

		double dl=0, dr=0; // delay of children
		double Cl=0, Cr=0;
		double dmax = 0;
		BtreeNode *lchd = rt->lchild, *rchd = rt->rchild;


		if (rt->getType()==LEAF) {
			// rt is a leaf
			dmax = 0;
			rt->dmax = 0;
			rt->dmin = 0;
			c = rt->getC();						
		} else if (lchd && rchd) {
			// rt is a merging point
			// traverse the children first
			dl = DME3D_simple(lchd, Cl, insert_buf); // determine Cl
			dr = DME3D_simple(rchd, Cr, insert_buf);

			if(insert_buf)
				 // the coordinate of rt is adjusted in this function, using Cl, Cr
				dmax = det_branch_between_nodes(rt, c, dl, Cl, dr, Cr);
			else
				dmax = det_branch_nobuf(rt, c, lchd, dl, Cl, rchd, dr, Cr); // no buffer is inserted
		} else if (rt->getType()!= VIA){
			// only TSV can have one child
			stringstream ss;
			ss << " Error: there is a node (" << rt->getType() << rt->getX() << "_"
				<< rt->getY() << ") in Btree with only one child ";
			if (!lchd)
				ss << rchd->getType() << rchd->getID() << "!";
			else
				ss << lchd->getType() << lchd->getID() << "!";
			throw(ss.str());
		} else {
			// a TSV with one child 
			double dmin;
			int x1, y1;
			BtreeNode * chd_temp=lchd;

			if(!lchd){
				chd_temp = rchd;
			}


			dmax = DME3D_simple(chd_temp, Cl, insert_buf);
			dmin = chd_temp->dmin;
			x1 = chd_temp->getX();
			y1 = chd_temp->getY();

			rt->setXY(x1, y1); // place TSV next to its child to gather buffers in the same tier

			if(insert_buf && Cl >= Buffer::cload_limit){
				// a buffer is needed right after TSV, uncomment if don't want buffer to be inserted
				Buffer *new_buf = new Buffer(num_buf++, x1,y1,chd_temp->getZ());
				insertNode(chd_temp, new_buf);
				Cl = new_buf->getC();
				new_buf->calc_path_delay();
			}


			c = rt->getC()+Cl;	// cload seen at the input of TSV
			rt->calc_path_delay(); // update delay
		}

		//cout << rt->getType() << " " << rt->dmax << " " << rt->dmin << endl;

		return dmax;
	}

	double Btree::det_branch_between_nodes(BtreeNode *rt, double& Ctotal, double dl, double Cl, 
		double dr, double Cr) {
			// determine the length of branches between lchd and rchd, dl is the largest delay from lchd to sinks, 
			// Cl is the cap seen at lchd, determine the locationi of rt
			double dmax = 0; // dmax is the largest delay from the rt to sinks of lchd and rchd
			double ddelta = dl - dr;
			double ddelta_wire = 0;
			double dl1, dr1;
			double r = Wire::R, c = Wire::C;
			int x, y, z;
			double Cl_new = 0, Cr_new = 0;
			BtreeNode *lchd = NULL, *rchd = NULL;


			if (!rt)
				throw(" Error: null node in Btree::det_branch_between_nodes!");
			else{
				lchd = rt->lchild;
				rchd = rt->rchild;
			}

			if (!lchd || !rchd)
				throw(" Error: null child in Btree::det_branch_between_nodes!");

			double l = *lchd - *rchd;
			double ll, lr;
			int xl = lchd->getX(), yl = lchd->getY();
			int xr = rchd->getX(), yr = rchd->getY();

			// ELMORE_DISTR* ll^2 * rc + ELMORE_LUMP * ll *r *Cl + ddelta = ELMORE_DISTR*lr^2*rc + ELMORE_LUMP*lr *r *Cr
			ll = (ELMORE_DISTR*r*c*pow(l,2) + ELMORE_LUMP*r*l*Cr - ddelta) / 
				(2*ELMORE_DISTR*r*c*l + ELMORE_LUMP* r * (Cl+Cr)); // to obtained 0-skew
			lr = l - ll;
			ddelta_wire = ELMORE_DISTR* Wire::R * Wire::C * pow(l,2); // largest delta d provided by the whole wire

			if (ll < 0){
				// detour or extra buffer is needed for lr, dl >> dr
				//cout << lchd->getType() << lchd->getID() << " ll : " << ll << endl
				double ddelta_buf = Buffer::Dbuf + 
					ELMORE_LUMP * (Buffer::Rout * Cr + Wire::R*l*Buffer::Cin); // d delta provide by a buffer

				if(dl >= dr+ddelta_buf+ddelta_wire){
					// insert a buffer right before right chd
					x = xr;
					y = yr;
					z = rchd->getZ();
					Buffer *new_buf = new Buffer(num_buf++, x, y, z);
					insertNode(rchd, new_buf); // insert new_buf before lchd
					new_buf->calc_path_delay();
					dr = new_buf->dmax;
					Cr = new_buf->getC();					
					dmax = det_branch_between_nodes(rt,Ctotal, dl, Cl, dr, Cr); // recalculate
					return dmax;
				} else{
					// tolerate a skew, buffer cannot help
					ll = 0;
					lr = l; // lr is simplified to l, detour wires are needed
					x = xl;
					y = yl;
				}
			} else if (lr < 0){
				// detour is needed for ll, dl << dr
				//cout << rchd->getType() << rchd->getID() << " lr : " << lr << endl;
				double ddelta_buf = Buffer::Dbuf + 
					ELMORE_LUMP * (Buffer::Rout * Cl + Wire::R*l*Buffer::Cin); // d delta provide by a buffer

				if(dr >= dl+ddelta_buf+ddelta_wire){
					// insert a buffer right before left chd
					x = xl;
					y = yl;
					z = lchd->getZ();
					Buffer *new_buf = new Buffer(num_buf++, x, y, z);
					insertNode(lchd, new_buf); // insert new_buf before lchd
					new_buf->calc_path_delay();
					dl = new_buf->dmax;
					Cl = new_buf->getC();					
					dmax = det_branch_between_nodes(rt,Ctotal, dl, Cl, dr, Cr); // recalculate
					return dmax;
				} else{
					// tolerate a skew, buffer cannot help
					ll = l;
					lr = 0; // ll is simplified to l, detour wires are needed
					x = xr;
					y = yr;					
				}
			} else {
				// rt(x,y) is in the middle of lchd and rchd
				x = (xr - xl)* ll/l + xl;
				y = (yr - yl)* ll/l + yl;
			}


			Cl_new = ll*Wire::C + Cl;
			Cr_new = lr*Wire::C + Cr;

			// check buffer load constraints
			if(Cl_new >= Buffer::cload_limit || Cr_new >= Buffer::cload_limit){
				if(Cl_new >= Buffer::cload_limit){
					// a buf is needed for the left branch
					int l_for_buf = max(Buffer::cload_limit - Cl, 0.0) / Wire::C; // max length that can be driven
					x = (xr-xl)*l_for_buf / l +xl;
					y = (yr-yl)*l_for_buf / l +yl;
					z = lchd->getZ();
					Buffer *new_buf1 = new Buffer(num_buf++, x, y, z);
					insertNode(lchd, new_buf1); // insert new_buf before lchd
					new_buf1->calc_path_delay();
					dl = new_buf1->dmax;
					Cl = new_buf1->getC();					
				}

				if(Cr_new >= Buffer::cload_limit){
					// a buf is needed for the right branch
					int l_for_buf = max(Buffer::cload_limit - Cr, 0.0) / Wire::C; // max length that can be driven
					x = (xl-xr)*l_for_buf / l +xr;
					y = (yl-yr)*l_for_buf / l +yr;
					z = rchd->getZ();
					Buffer *new_buf2 = new Buffer(num_buf++, x, y, z);
					insertNode(rchd, new_buf2); // insert new_buf before lchd
					new_buf2->calc_path_delay();
					dr = new_buf2->dmax;
					Cr = new_buf2->getC();
				}

				dmax = det_branch_between_nodes(rt,Ctotal, dl, Cl, dr, Cr); // recalculate
				return dmax;
			}
			//else if(Cl_new + Cr_new > Buffer::cload_limit){
			//	// the sum exceeds but chd not
			//} 
			else {
				// buffer load constraint is satisfied, setting rt locations and delay
				rt->setXYZ(x, y, rt->getZ()); // z doesn't change

				dl1 = ELMORE_DISTR * pow(ll, 2) * r * c + ELMORE_LUMP * r * ll * Cl;
				dr1 = ELMORE_DISTR * pow(lr, 2) * r * c + ELMORE_LUMP * r * lr * Cr;

				Ctotal = l * c + Cl + Cr;
				dmax = max(dl1+dl, dr1+dr) + ELMORE_LUMP * rt->R * Ctotal + rt->getD();
				if (rt->getType() != BUF)
					Ctotal = Ctotal + rt->getC();
				else
					Ctotal = rt->getC();

				rt->dmax = dmax;
				rt->dmin = min(lchd->dmin + dl1, rchd->dmin + dr1);
			}

			return dmax;
	}

	void Btree::CTS() {
		// clock tree synthesis for 2D
		root = MMM(sinks);
	}

	void Btree::setRoot(BtreeNode *elem) {
		// set the content of the root inside "first"
		root = elem;
		setParent(root, NULL);
	}

	//
	bool Btree::isEmpty() {
		if (nsink == 0)
			return true;
		else
			return false;
	}

	void Btree::setParent(BtreeNode *elem, BtreeNode *par) {
		// set the parent of elem to par
		if (elem)
			// elem is not null
			elem->parent = par;
	}

	BtreeNode* Btree::centerofset(const vector<BtreeNode*>& nodes) {
		// return the among the nodes
		int n = nodes.size();
		int x = 0;
		int y = 0;
		int z = 0;
		BtreeNode* cnode = NULL;

		if (n == 1)
			return nodes[0];
		else if (n == 0)
			return NULL;

		for (int i = 0; i <= n - 1; i++) {
			x = x + nodes[i]->x;
			y = y + nodes[i]->y;
		}

		x = x / n;
		y = y / n;
		z = -1; // to be replaced
		cnode = new BtreeNode(-2, CON, x, y, z); // an internal node, -2 for debugging

		return cnode;

	}

	vector<Cxuh::BtreeNode*> Btree::subSet(const vector<BtreeNode*>& nodes,
		BtreeNode* center, char direction) {
			// divide a set into two sets
			vector<Cxuh::BtreeNode*> snodes;
			int x0 = center->getX();
			int y0 = center->getY();
			int i;
			int n = nodes.size();

			if (n <= 1)
				// if nodes.size == 1, return a empty subset
				return snodes;

			for (i = 0; i <= n - 1; i++) {
				switch (direction) {
				case 'l':
					if (nodes[i]->getX() <= x0)
						snodes.push_back(nodes[i]);
					break;
				case 'r':
					if (nodes[i]->getX() > x0)
						snodes.push_back(nodes[i]);
					break;
				case 'b':
					if (nodes[i]->getY() <= y0)
						snodes.push_back(nodes[i]);
					break;
				case 't':
					if (nodes[i]->getY() > y0)
						snodes.push_back(nodes[i]);
					break;
				default:
					cerr << "Wrong MMM subset!" << endl;
					break;
				}
			}

			return snodes;

	}

	int Btree::d2subcenter(const vector<BtreeNode*>& nodes, char direction) {
		// delay from center of nodes to the center of its subset, in the specified direction
		int d = 0;
		BtreeNode* cnode = centerofset(nodes);
		BtreeNode* scnode; // center node of subset
		vector<BtreeNode*> snodes;

		switch (direction) {
		case 'x':
			// distance to left subset, equal to right subset
			snodes = subSet(nodes, cnode, 'l');
			break;
		case 'y':
			snodes = subSet(nodes, cnode, 't');
			break;
		default:
			cerr << "Wrong direction of d2subcenter()!" << endl;
			break;
		}

		scnode = centerofset(snodes); // if snodes is empty, return NULL

		if (!cnode && !scnode)
			d = (*cnode) - (*scnode);
		else
			d = numeric_limits<int>::max(); // nodes or snodes is empty, return max distance

		return d;
	}

	BtreeNode* Btree::MMM(vector<BtreeNode*>& nodes) {
		// top-down MMM
		BtreeNode* sroot = NULL;
		double dmax = numeric_limits<double>::max();
		double d1, d2; // delay
		double s1 = dmax, s2 = dmax; // initialize skew to maximum
		vector<BtreeNode*> Sl, Sr, St, Sb;

		if (nodes.size() <= 1) {
			// only one node
			sroot = nodes[0];
			return sroot;
		}

		sroot = centerofset(nodes);
		if (sroot->getType() != LEAF) {
			//a merge point is generated
			sroot->id = num_mp;
			num_mp += 1; // a new merge point is generated
		}

		Sl = subSet(nodes, sroot, 'l'); // left subset
		Sr = subSet(nodes, sroot, 'r'); // right subset
		St = subSet(nodes, sroot, 't'); // top subset
		Sb = subSet(nodes, sroot, 'b'); // bottom subset

		// delay equalization look-ahead
		if (Sl.size() && Sr.size()) {
			// Sl and Sr are not empty
			d1 = d2subcenter(Sl, 'y');
			d2 = d2subcenter(Sr, 'y');
			s1 = abs(d1 - d2);
		} else if (St.size() && Sb.size()) {
			// St and Sb are not empty
			d1 = d2subcenter(St, 'x');
			d2 = d2subcenter(Sb, 'x');
			s2 = abs(d1 - d2);
		} else
			cerr << "Cannot divide nodes!" << endl;

		if (s1 <= s2) {
			// use horizontal division for current set
			// vertical division for subsets
			sroot->lchild = MMM(Sl, 'y');
			sroot->rchild = MMM(Sr, 'y');
		} else {
			// use vertical division for current set
			// horizontal division for subsets
			sroot->lchild = MMM(St, 'x');
			sroot->rchild = MMM(Sb, 'x');
		}

		if (sroot->lchild)
			setParent(sroot->lchild, sroot);
		if (sroot->rchild)
			setParent(sroot->rchild, sroot);

		return sroot;
	}

	BtreeNode* Btree::MMM(vector<Cxuh::BtreeNode*>& nodes, char direction) {
		// top-down MMM with specified direction
		BtreeNode* sroot = NULL;
		vector<BtreeNode*> Sl, Sr, St, Sb;

		if (nodes.size() == 1) {
			// only one node
			sroot = nodes[0];
			return sroot;
		} else if (nodes.size() == 0) {
			return NULL;
		}

		sroot = centerofset(nodes);
		if (sroot->getType() != LEAF) {
			//a merge point is generated
			sroot->id = num_mp;
			num_mp += 1; // a new merge point is generated
		}

		if (direction == 'x') {
			// use horizontal division for current set
			// determine by delay one level look-ahead for subsets
			Sl = subSet(nodes, sroot, 'l'); // left subset
			Sr = subSet(nodes, sroot, 'r'); // right subset
			if (Sl.size() && Sr.size()) {
				sroot->lchild = MMM(Sl);
				sroot->rchild = MMM(Sr);
			} else
				cerr << "Try to divide a set in wrong direction!" << endl;
		} else {
			// use vertical division for current set
			// determine by delay one level look-ahead for subsets
			St = subSet(nodes, sroot, 't'); // top subset
			Sb = subSet(nodes, sroot, 'b'); // bottom subset
			if (St.size() && Sb.size()) {
				sroot->lchild = MMM(St);
				sroot->rchild = MMM(Sb);
			} else
				cerr << "Try to divide a set in wrong direction!" << endl;
		}

		if (sroot->lchild)
			setParent(sroot->lchild, sroot);
		if (sroot->rchild)
			setParent(sroot->rchild, sroot);

		return sroot;
	}

	BtreeNode* Btree::MMM3D(vector<Cxuh::BtreeNode*>& nodes, int TSVbound,
		char direction) {
			// top-down MMM considering Z-cut and TSV bound (min #TSV within one tier)
			BtreeNode* sroot = NULL, *sroot1 = NULL, *sroot2 = NULL;
			vector<Cxuh::BtreeNode*> nodes1, nodes2; // subsets of nodes
			int TSVbound1 = 1, TSVbound2 = 1;
			double TSVratio;
			char direction2 = direction;
			int zmin, zmax;

			// check whether nodes are leaves
			if (nodes.empty()) {
				throw( " Error: MMM3D for EMPTY nodes!");
			} else if (nodes.size() == 1)
				return nodes[0];

			// check z of nodes
			findzMinMax(nodes, zmin, zmax);

			if (TSVbound == 0)
				throw(" Error: TSV bound for MMM3D is 0!");
			else if (TSVbound == 1 && zmin != zmax) {
				//only one tsv can be used and the nodes are in different tiers, Z cut for 3D nodes
				sroot = Zcut(nodes, nodes1, nodes2);
				TSVbound1 = 1;
				TSVbound2 = 1;
			} else {
				// XY cut, 2D MMM
				sroot = XYcut(nodes, nodes1, nodes2, TSVratio, direction); // ratio = TSVbound1/TSVbound2
				TSVbound1 = max(Txuh::round(TSVbound * (TSVratio / (1 + TSVratio))), 1); // must be larger or equal to 1
				TSVbound2 = TSVbound - TSVbound1;
				if (TSVbound2 == 0) {
					TSVbound2 = 1;
					TSVbound1 = max(TSVbound - TSVbound2, 1); // TSVbound2 has been reset, TSVbound1 needs to be re-calculated again
				}

				if (direction == 'x')
					direction2 = 'y';
				else
					direction2 = 'x';
			}

			sroot1 = MMM3D(nodes1, TSVbound1, direction2);
			sroot2 = MMM3D(nodes2, TSVbound2, direction2);

			// generate TSVs if needed between sroot and sroot1
			sroot->lchild = sroot1;
			sroot->rchild = sroot2;
			setParent(sroot1, sroot);
			setParent(sroot2, sroot);
			genTSV(sroot, sroot1);
			genTSV(sroot, sroot2);

			return sroot;
	}

	BtreeNode* Btree::XYcut(const vector<Cxuh::BtreeNode*>& nodes,
		vector<Cxuh::BtreeNode*>& nodes1, vector<Cxuh::BtreeNode*>& nodes2,
		double& TSVratio, char direction) {
			// XY cut for 3D nodes, return sroot, nodes1, nodes2, ratio = TSVbound1/TSVbound2
			BtreeNode* sroot = NULL;
			int zmin = numeric_limits<int>::max(), zmax = 0;
			int nmin1 = 0, nmin2 = 0;

			if (nodes.size() == 1) {
				// only one node
				sroot = nodes[0];
				return sroot;
			} else if (nodes.size() == 0) {
				// empty nodes
				return NULL;
			}

			// generate a sroot and mark it as merge point if true
			sroot = centerofset(nodes);
			if (sroot->getType() != LEAF) {
				//a merge point is generated
				sroot->id = num_mp;
				num_mp += 1; // a new merge point is generated
			}

			//determine z of sroot
			findzMinMax(nodes, zmin, zmax);
			if (zmin > srcz)
				sroot->z = zmin;
			else if (zmax < srcz)
				sroot->z = zmax;
			else
				sroot->z = srcz;

			if (direction == 'x') {
				// use horizontal division for current set
				nodes1 = subSet(nodes, sroot, 'l'); // left subset
				nodes2 = subSet(nodes, sroot, 'r'); // right subset
				if (nodes1.size() == 0 || nodes2.size() == 0) {
					// x of all nodes are the same, y-cut should be tried
					nodes1 = subSet(nodes, sroot, 't'); // top subset
					nodes2 = subSet(nodes, sroot, 'b'); // bottom subset
				}
			} else {
				// use vertical division for current set
				nodes1 = subSet(nodes, sroot, 't'); // top subset
				nodes2 = subSet(nodes, sroot, 'b'); // bottom subset
				if (nodes1.size() == 0 || nodes2.size() == 0) {
					// y of all nodes are the same, x-cut should be tried
					nodes1 = subSet(nodes, sroot, 'l'); // top subset
					nodes2 = subSet(nodes, sroot, 'r'); // bottom subset
				}
			}

			// determine TSV bound ratio
			nmin1 = minN(nodes1);
			nmin2 = minN(nodes2);

			if (nmin1 == 0 || nmin2 == 0) {
				cerr << "Min # of nodes is 0 for a set of divided nodes!" << endl;
				return NULL;
			}
			TSVratio = double(nmin1) / double(nmin2);

			return sroot;
	}

	int Btree::minN(const vector<Cxuh::BtreeNode*>& nodes) {
		// min nsink within the same tier in nodes
		int n = numeric_limits<int>::max(), nodesize = nodes.size(), z;
		map<int, int> ns; // nsinks in different tiers
		map<int, int>::iterator it;

		// determine n for each tier
		for (int i = 0; i <= nodesize - 1; i++) {
			z = nodes[i]->getZ();
			if (ns.count(z)) {
				// similar sinks have been found
				ns[z] = ns[z] + 1;
			} else
				// insert new tier
				ns.insert(pair<int, int>(z, 1));
		}

		// find the min n
		for (it = ns.begin(); it != ns.end(); it++) {
			if ((*it).second < n)
				n = (*it).second;
		}

		return n;

	}

	BtreeNode* Btree::Zcut(const vector<Cxuh::BtreeNode*>& nodes,
		vector<Cxuh::BtreeNode*>& nodes1, vector<Cxuh::BtreeNode*>& nodes2) {
			// Z cut for 3D nodes, return sroot, nodes1, nodes2
			BtreeNode* sroot = NULL;
			int n = nodes.size();
			int zmin = numeric_limits<int>::max(), zmax = 0;
			int zstd = 0;

			if (n <= 1) {
				nodes1.clear();
				nodes2.clear();
				return sroot; // cannot be divided
			}

			// generate a sroot and mark it as merge point if true
			sroot = centerofset(nodes);
			if (sroot->getType() != LEAF) {
				//a merge point is generated
				sroot->id = num_mp;
				num_mp += 1; // a new merge point is generated
			}

			findzMinMax(nodes, zmin, zmax);

			// determine which z should be used to divide nodes and sroot should be
			if (zmin > srcz)
				zstd = zmin;
			else if (zmax < srcz)
				zstd = zmax;
			else
				zstd = srcz;

			sroot->z = zstd;

			// divides nodes
			for (int i = 0; i <= n - 1; i++) {
				int iz = nodes[i]->z;
				if (iz == zstd)
					nodes1.push_back(nodes[i]);
				else
					nodes2.push_back(nodes[i]); // assuming there are already nodes in srcz tier
			}

			// there is no node at zstd, divide in another way
			if (nodes1.empty()) {
				nodes2.clear();
				for (int i = 0; i <= n - 1; i++) {
					int iz = nodes[i]->z;
					if (iz < zstd)
						nodes1.push_back(nodes[i]);
					else
						nodes2.push_back(nodes[i]); // assuming there are already nodes in srcz tier
				}
			}

			return sroot;
	}

	void Btree::findzMinMax(const vector<Cxuh::BtreeNode*>& nodes, int& min1,
		int& max1) {
			// find the min&max z of a set
			int n = nodes.size();

			if (!n) {
				cerr << "Try to find z min&max for an empty node set!" << endl;
				return;
			}

			// initialization
			min1 = numeric_limits<int>::max();
			max1 = 0;

			for (int i = 0; i <= n - 1; i++) {
				int iz = nodes[i]->z;
				if (iz < min1)
					min1 = iz;
				if (iz > max1)
					max1 = iz;
			}

	}

	void Btree::genTSV(BtreeNode* sroot, BtreeNode* sroot1) {
		// generate a TSV between sroot and sroot1, child = 'r' or 'l'
		if (sroot == NULL || sroot1 == NULL)
			return;

		TSV* newtsv = NULL;
		//int x0 = sroot->getX();
		int x1 = sroot1->getX();
		//int y0 = sroot->getY();
		int y1 = sroot1->getY();
		int z0 = sroot->getZ();
		int z1 = sroot1->getZ();

		if (z0 != z1) {
			// TSVs are needed between sroot and sroot1
			newtsv = new TSV(num_tsv, VIA, x1, y1, z0, z1); // tsv is placed close to the sinks

			num_tsv += 1;
			if (sroot1->getType() == CON)
				// sroot1 is a merge point, can be deleted
				replaceNode(sroot1, newtsv);
			else
				// sroot1 is a buf, tsv, or leaf, can't be deleted
				insertNode(sroot1, newtsv);
		}

	}

	void Btree::replaceNode(BtreeNode* oldnd, BtreeNode* newnd) {
		// replace oldnd with newnd
		BtreeNode* temp = NULL;

		if (!oldnd || !newnd)
			return;

		temp = oldnd->parent;
		setParent(newnd, temp);
		if (temp != NULL) {
			if (temp->lchild == oldnd)
				temp->lchild = newnd;
			else
				temp->rchild = newnd;
		}

		newnd->lchild = oldnd->lchild;
		newnd->rchild = oldnd->rchild;

		setParent(newnd->lchild, newnd);
		setParent(newnd->rchild, newnd);

		delete oldnd;
	}

	void Btree::insertNode(BtreeNode* oldnd, BtreeNode* newnd) {
		// insert newnd before oldnd
		BtreeNode* temp = NULL;

		if (!oldnd || !newnd)
			return;

		temp = oldnd->parent;
		setParent(newnd, temp);
		// connect parent with newnd
		if (temp != NULL) {
			if (temp->lchild == oldnd)
				temp->lchild = newnd;
			else
				temp->rchild = newnd;
		}

		// insert oldnd to the empty child of newnd
		if (!newnd->lchild)
			newnd->lchild = oldnd;
		else if (!newnd->rchild)
			newnd->rchild = oldnd;
		else {
			string errorstr = "Try to insert a subtree between two nodes!";
			throw(errorstr);
		}
		setParent(oldnd, newnd);
	}

	void Btree::buffering(BtreeNode* rt, double Cprev) {
		// top-down slew-aware buffer insertion for the tree rooted from rt, no downstream buffer!
		// constraint is Buffer::cload_limit
		// the cap seen by the parent of rt is Cprev

		if (!rt)
			return;

		if(Cprev > Buffer::cload_limit || Cprev < 0){
			string errstr = " Error: while buffering, the load of previous stage is out of range!";
			throw(errstr);
		}

		if (Buffer::cload_limit <= Buffer::Cin)
			throw("Buffer::cload_limit is smaller than Buffer::Cin!");


		double c1 = 0, c2 = 0;	// cap between rt and its children
		BtreeNode *chdl = rt->lchild, *chdr = rt->rchild;
		double cload = 0; // load of the parent of rt
		int l; // distance between buffers and rt
		int l1, l2; // distance between rt and its children

		if(chdl && chdr){
			l1 = *rt - *chdl;
			l2 = *rt - *chdr;
			c1 = Wire::C*l1 + chdl->getC();
			c2 = Wire::C*l2 + chdr->getC();
			cload = Cprev + c1 + c2;
			if(cload >= Buffer::cload_limit){
				// a buffer is needed
				l = (Buffer::cload_limit - Cprev)/Wire::C/2;
				l1 = min(l1, l);
				l2 = min(l2, l);
			}
		} else if (chdl){
			c1 = Wire::C*(*rt - *chdl) + chdl->getC();
		} else if (chdr){
			c2 = Wire::C*(*rt - *chdr) + chdr->getC();
		} else
			return;	// no children, no buffer is needed
	}

	void Btree::buffering(BtreeNode* rt) {
		// slew-aware buffer insertion for the tree rooted from rt

		if (!rt)
			return;

		if (Buffer::cload_limit <= Buffer::Cin)
			throw("Buffer::cload_limit is smaller than Buffer::Cin!");

		double cload;
		int l, l1, x1, y1, z1;
		Buffer* newbuf = NULL;
		BtreeNode* par = rt->parent;

		buffering(rt->lchild);
		buffering(rt->rchild);

		cload = getInload(rt);

		if(cload<0)
			throw("GetInload() returns negative C!");

		l = max(Buffer::cload_limit - cload, 0.0) / Wire::C; // max length that can be driven

		if (!par) {
			// rt is the root of the tree, a buffer is insert before it by default
			x1 = rt->x;
			y1 = rt->y;
			z1 = rt->z;
			newbuf = new Buffer(num_buf, x1, y1, z1, Buffer::Rout, Buffer::Cin, Buffer::Dbuf, Buffer::cload_limit);
			replaceNode(rt, newbuf);
			root = newbuf;
		} else {
			l1 = *rt - *par; // distance between rt and its parent
			while (l1 > l / 2) {
				// a buffer is needed
				if (l1 > l) {
					// insert a buffer between par and rt
					double dx = fabs(float(par->x - rt->x));
					double dy = fabs(float(par->y - rt->y));

					if (l <= dx && dx > 0) {
						// horizontally extending rt
						x1 = double(l) / dx * (par->x - rt->x) + rt->x;
						y1 = rt->y;
					} else if (dx ==0 && dy == 0) {
						throw("Inserting buffer into two nodes at the same location!");
					} else {
						// first horizontal then vertical
						x1 = par->x;
						y1 = (l - dx) / dy * (par->y - rt->y) + rt->y;
					}

					z1 = rt->z;
					l1 = l1 - l;
				} else {
					// a buffer right after the parent
					x1 = par->x;
					y1 = par->y;
					z1 = par->z;
					l1 = 0;
				}
				if(x1 < 0 || y1 < 0 || z1 < 0){
					stringstream errss;
					errss<< "The x,y,z of a buffer is negative! " << x1 << "," << y1 << "," << z1;
					throw(errss.str());
				}

				newbuf = new Buffer(num_buf, x1, y1, z1, Buffer::Rout, Buffer::Cin, Buffer::Dbuf, Buffer::cload_limit);
				insertNode(rt, newbuf); // insert new buf before rt
				rt = newbuf; // check whether a buf is needed before newbuf
				l = (Buffer::cload_limit - Buffer::Cin) / Wire::C;
			}
		}

	}

	double Btree::getInload(BtreeNode* nd) {
		// get the cap load seen at the input of nd, including its children
		return nd->get_in_cload();
	}

	void Btree::printPaths(std::ostream& out) {
		// print all paths
		unsigned int i = 0;
		int l = 0;
		BtreeNode* tempnd, *par;
		string t;

		for (i = 0; i <= sinks.size() - 1; i++) {
			tempnd = sinks[i];
			l = 0;
			while (tempnd) {
				par = tempnd->parent;
				out << tempnd->id << "(";
				int x1 = tempnd->x;
				int y1 = tempnd->y;
				int z1 = tempnd->z;

				switch (tempnd->ntype) {
				case CON:
					t = "Merge";
					break;
				case LEAF:
					t = "Sink";
					break;
				case VIA:
					t = "TSV";
					break;
				case BUF:
					t = "BUF";
					break;
				default:
					t = "Other";
					break;
				}
				if(x1 < 0 || y1 < 0 || z1 < 0){
					stringstream errss;
					errss<< "The x,y,z of a buffer is negative! " << x1 << "," << y1 << "," << z1;
					throw(errss.str());
				}

				out << t << ", " << tempnd->z << "): " << getInload(tempnd);
				if (par != NULL) {
					out << " --" << *tempnd - *par << "--> "; // d between nd and par
					l = l + (*tempnd - *par);
				}
				tempnd = par;
			}
			out << ". Length: " << l;
			out << endl;
		}
	}

	void Btree::findXYMinMax(int& xmin, int& xmax, int& ymin, int& ymax) {
		// get the boundary of sinks
		unsigned int i;
		BtreeNode* temp;

		xmin = numeric_limits<int>::max();
		xmax = 0;
		ymin = xmin;
		ymax = 0;

		for (i = 0; i <= sinks.size() - 1; i++) {
			temp = sinks[i];
			if (temp->x > xmax)
				xmax = temp->x;
			if (temp->x < xmin)
				xmin = temp->x;
			if (temp->y > ymax)
				ymax = temp->y;
			if (temp->y < ymin)
				ymin = temp->y;
		}
	}

	void Btree::display4sim(ostream& out, const int id1) {
		// display shared and independent segments for sim
		int l1;
		BtreeNode *temp1=NULL, *par1=NULL;
		string t;

		temp1 = sinks[id1];

		l1 = 0;
		while (temp1) {
			par1 = temp1->parent;
			out << temp1->id << "(";

			out << temp1->z << "): " << l1 << " -- ";
			if (par1 != NULL) {
				l1 = l1 + (*temp1 - *par1);
			}
			temp1 = par1;
		}
		out << endl;
	}

	void Btree::DME3D(){
		// DME 3D including buffer insertion
		MergingSegment* ms = DME3D_1(root);	// bottom-up merging
		root = DME3D_2(ms);	// top-down buf insertion
	}

	MergingSegment* Btree::DME3D_1(BtreeNode* nd){
		// DME 3D 1st phase, bottom-up
		MergingSegment* ms = NULL, *lms = NULL, *rms = NULL;
		BtreeNode* lchd, *rchd;

		if (!nd)
			throw(" Error: DME3D_1 on a NULL!");
		else{
			lchd = nd->lchild;
			rchd = nd->rchild;
		}

		if (!lchd && !rchd){
			// nd is a leaf
			ms = new MergingSegment(nd);
		}else if (!lchd || !rchd) {
			throw(" Btree::DME3D_1(): There is a node with only one child!");
		}else{
			// with two leaves
			lms = DME3D_1(lchd);
			rms = DME3D_1(rchd);
			ms = new MergingSegment(nd, lms, rms);
		}

		return ms;
	}

	BtreeNode* Btree::DME3D_2(MergingSegment* ms){
		// DME 3D 2nd phase, top-down
		BtreeNode* nd = NULL;

		return nd;
	}

	void Btree::update(){
		// update the count and label of all nodes
		num_buf = 0;
		num_mp = 0;
		num_tsv = 0;
		pd_min = numeric_limits<double>::max();
		pd_max = numeric_limits<double>::min();

		update(root);
	}

	void Btree::update(BtreeNode* rt){
		// update the count and label in the subtree
		if(!rt)
			return;	// empty node

		stringstream ss;

		if(rt->ntype == BUF){
			rt->id = num_buf;
			num_buf++;
		} else if (rt->ntype == VIA){
			rt->id = num_tsv;
			num_tsv++;
		} else if (rt->ntype == CON){
			rt->id = num_mp;
			num_mp++;
		}

		ss << rt->ntype << rt->id;
		rt->name = ss.str();

		if (rt->lchild)
			update(rt->lchild);
		if (rt->rchild)
			update(rt->rchild);

		rt->calc_path_delay(); // update dmax and dmin
	}

	double Btree::get_pathd_min() const{
		// get min delay
		if (!root)
			return numeric_limits<double>::max();
		else
			return root->dmin;
	}

	double Btree::get_pathd_max() const{
		// get max delay
		if (!root)
			return numeric_limits<double>::min();
		else
			return root->dmax;
	}

	double Btree::det_branch_nobuf(BtreeNode *rt, double& Ctotal,
		BtreeNode *lchd, double dl, double Cl, BtreeNode *rchd, double dr,
		double Cr) {
			// determine the length of branches between lchd and rchd, dl is the largest delay from lchd to sinks, 
			// Cl is the cap seen at lchd, determine the locationi of rt
			// no buffer is inserted
			double dmax = 0; // dmax is the largest delay from the rt to sinks of lchd and rchd
			double ddelta = dl - dr;
			double dl1, dr1;
			double r = Wire::R, c = Wire::C;
			int x, y;
			double Cl_new = 0, Cr_new = 0;

			if (!rt || !lchd || !rchd)
				throw(" Error: null node in Btree::det_branch_between_nodes!");

			double l = *lchd - *rchd;
			double ll, lr;
			int xl = lchd->getX(), yl = lchd->getY();
			int xr = rchd->getX(), yr = rchd->getY();

			// ELMORE_DISTR* ll^2 * rc + ELMORE_LUMP * ll *r *Cl + ddelta = ELMORE_DISTR*lr^2*rc + ELMORE_LUMP*lr *r *Cr
			ll = (ELMORE_DISTR*r*c*pow(l,2) + ELMORE_LUMP*r*l*Cr - ddelta) / 
				(2*ELMORE_DISTR*r*c*l + ELMORE_LUMP* r * (Cl+Cr)); // to obtained 0-skew
			lr = l - ll;

			if (ll < 0){
				// detour is needed for lr
				//cout << lchd->getType() << lchd->getID() << " ll : " << ll << endl;
				ll = 0;
				lr = l; // lr is simplified to l, detour wires are needed
				x = xl;
				y = yl;
			} else if (lr < 0){
				// detour is needed for ll
				//cout << rchd->getType() << rchd->getID() << " lr : " << lr << endl;
				lr = 0;
				ll = l;
				x = xr;
				y = yr;
			} else {
				// rt(x,y) is in the middle of lchd and rchd
				x = (xr - xl)* ll/l + xl;
				y = (yr - yl)* ll/l + yl;
			}

			Cl_new = ll*Wire::C + Cl;
			Cr_new = lr*Wire::C + Cr;

			
			// buffer load constraint is satisfied, setting rt locations and delay
			rt->setXYZ(x, y, rt->getZ()); // z doesn't change

			dl1 = ELMORE_DISTR * pow(ll, 2) * r * c + ELMORE_LUMP * r * ll * Cl;
			dr1 = ELMORE_DISTR * pow(lr, 2) * r * c + ELMORE_LUMP * r * lr * Cr;

			Ctotal = l * c + Cl + Cr;
			dmax = max(dl+dl1, dr+dr1) + ELMORE_LUMP * rt->R * Ctotal + rt->getD();
			if (rt->getType() != BUF)
				Ctotal = Ctotal + rt->getC();
			else
				Ctotal = rt->getC();

			rt->dmax = dmax;
			rt->dmin = min(lchd->dmin + dl1,
				rchd->dmin + dr1);
			
			return dmax;
	}

}
